# !diagnostics off
library(readr)
library(magrittr)
library(sandwich)
library(xtable)
library(ggplot2)
library(lmtest)
library(haven)
library(fixest)
library(Hmisc)
library(stargazer)
library(vtable)
library(conleyreg)
library(plm)
library(dplyr)
library(rdrobust)
library(descr)
library(ggthemes)
library(interplot)
library(mediation)
library(readxl)
library(broom)
library(tidyverse)
library(hrbrthemes)
library(logistf) 
library(gridExtra)
library(dotwhisker)
library(sf)
library(spdep)
library(spatialreg)
library(splm)
library(geosphere)
library(broom.mixed)
library(ggpubr)
library(jtools)
library(tigris)
library(lfe)
library(devtools)
library(survival)
library(ggsurvfit)
library(ggfortify)


##set working directory
  # setwd(getwd())

### must run zero policy options twice; otherwise won't work 
spdep::set.ZeroPolicyOption(TRUE)
spdep::set.ZeroPolicyOption(TRUE)
spatialreg::set.ZeroPolicyOption(TRUE)
spatialreg::set.ZeroPolicyOption(TRUE)
options(scipen= 9)


###### Call in data #################
base_exp <- read_csv("main-gridyear.csv")
base <- read_csv("main-grid.csv")
base_exp %<>% arrange(gridid, year)


### reorder variables for analysis (panel)
base_exp =base_exp  %>% mutate( 
  capord= factor(capord, ordered = F, levels = c( 'Low','Med',  'High')), 
  capord_orig= factor(capord_orig, ordered = F, levels = c('Low', 'Med', 'High')),
  capord_four= factor(capord_four, ordered = F, levels = c('None','Low', 'Med', 'High')),
  leadcapord= factor(leadcapord, ordered = F, levels = c('Low', 'Med', 'High')),
  leadcapord_noNA= factor(leadcapord_noNA, ordered = F, levels = c('Low', 'Med', 'High'))
)
uoabord_exp <- base_exp %>% filter(uoabord ==1 )


### reorder variables for analysis (non-panel)
base  %<>%  mutate_at(names(base)[grepl('capord$' , names(base))],
                     ~factor(., ordered = F, levels = c( 'Low','Med',  'High'))                     )
uoabord = base %>% filter(uoabord==1)               


###### Make into panel data ########
## All gridcells:
 
base_exp <- ungroup(base_exp)
pbase= pdata.frame(base_exp, index = c('gridid', 'year'))

pbase<- pbase %>% dplyr::select('gridid', 'year', everything()) 
pbase$realyear <- pbase$year %>% as.numeric()+1829
pbase<- pbase %>% arrange(gridid, year)

## Highrisk gridcells:


puoab= pdata.frame(uoabord_exp, index = c('gridid', 'year'))
puoab<- puoab %>% dplyr::select('gridid', 'year', everything()) 
puoab$realyear <- puoab$year %>% as.numeric()+1829
puoab<- puoab %>% arrange(gridid, year) 

 
##### Check data ###

## Check balance: should be true
 is.pbalanced(pbase)
 is.pbalanced(puoab)

## check basic summary stats 
base_exp$capord%>% table(useNA= 'ifany') 
#    Low     Med    High 
# 1200430   59288   37850 
uoabord_exp$capord %>% table(useNA= 'ifany') 
# Low    Med   High 
# 648738  36351  26819 
base_exp$disputed%>% table(useNA= 'ifany') 
# 0       1 
# 1176034  121534 
uoabord_exp$disputed%>% table(useNA= 'ifany') 
# 0       1 
# 601287 110621 
base_exp$init%>% table(useNA= 'ifany')  
#       0       1    <NA> 
# 1176034    2122  119412 
uoabord_exp$init%>% table(useNA= 'ifany')  
#      0      1   <NA> 
# 601287   1925 108696 

## should be True 
all.equal(base_exp$disputed, pbase$disputed %>% as.numeric)
all.equal(uoabord_exp$disputed, puoab$disputed %>% as.numeric)




#-------------------------------------------------------------------------
######### Make neighboring matrices for spatial analysis ###########
#-------------------------------------------------------------------------

## call in spatial information 
spatial <- read_sf('rexport.shp', stringsAsFactors = F)
# leave only gridid and geometry (get rid of unnecessary varialbes)
geodat <- spatial %>%  dplyr::select(OBJECTID_2,geometry)

# Combine spatial and non-spatial data
visual <- tigris::geo_join(
  spatial_data = geodat,
  data_frame = base,
  by_sp = "OBJECTID_2",
  by_df = "gridid",
  how = "inner"
)

####  King's matrix for All Gridcells 
# style W means row-standardized 
nb<- poly2nb(visual, queen = F)
spatial_w <- nb2listw(nb, style = 'W', zero.policy = T)  
summary(spatial_w, zero.policy=T)
# Number of nonzero links: 29470 
# Percentage nonzero weights: 0.05178176 
# Average number of links: 3.906416 

###### subset to High risk areas 
uoabord_vis <- visual %>% filter(uoabord==1)
uoabord_vis$gridid = uoabord_vis$OBJECTID_2

###  King's matrix for Highrisk Gridcells 
nb<- poly2nb(uoabord_vis, queen= F)
spatial_w_uoab <- nb2listw(nb, style = 'W', zero.policy = T)  




 
##########################################################
############## Main Text Figures and Tables ##############
##########################################################
 

#-------------------------------------------------------
##  Figure 4 ######
#-------------------------------------------------------


visual_exp <- tigris::geo_join(
  spatial_data = geodat,
  data_frame = base_exp,
  by_sp = "OBJECTID_2",
  by_df = "gridid",
  how = "inner"
)

formap<- visual_exp %>% filter(year==2001) %>%  
  ggplot(aes(fill = as.factor(forfish))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2' ,  'chocolate4') , name= '113-115 \nForestry & Fish')

cropmap<- visual_exp %>% filter(year==2001) %>%  
  ggplot(aes(fill = as.factor(cropland))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2' ,  'chocolate4') , 
                    name= '111-112 \nCrop & Animal \nProduction')
oilmap<- visual_exp %>% filter(year==2001) %>%  
  ggplot(aes(fill = as.factor(oilpresent))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2' ,  'chocolate4') , name= '211 \nOil & Gas \nExtraction')
minemap<- visual_exp %>% filter(year==2001) %>%  
  ggplot(aes(fill = as.factor(minpres_all))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2' ,  'chocolate4') , name= '212 \nMining \n(Except Oil & \nGas)')
a<- gridExtra::grid.arrange( cropmap, formap,   oilmap,minemap,
                             ncol=2, nrow=2)
 
# ggsave(plot=a, 'resources_color.png', width = 12, height = 8)

 

#-------------------------------------------------------
###### Figure 5 ######
#-------------------------------------------------------


a<- ggplot2::ggplot(visual_exp %>% filter(year==1830), aes(fill = as.factor(capord))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2',  'orangered1',  'chocolate4'), 
                    name= 'Resource  Capital Intensity (Ordinal)')
b<- ggplot2::ggplot(visual_exp %>% filter(year==1880), aes(fill = as.factor(capord))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2',  'orangered1',  'chocolate4'), 
                    name= 'Resource  Capital Intensity  (Ordinal)')  

d<-  ggplot2::ggplot(visual_exp %>% filter(year==1920), aes(fill = as.factor(capord))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2',  'orangered1',  'chocolate4'), 
                    labels =c('Low', "Med (Minerals)", "High (Oil)"),
                    name= 'Resource \nCapital Intensity \n(Potential Benefit \nConcentration)') +
  theme(legend.title = element_text(size = 12, face = 'bold'),
        legend.text = element_text(size = 14))   

f=  ggplot2::ggplot(visual_exp %>% filter(year==1950), aes(fill = as.factor(capord))) + 
  geom_sf(color= '#FFC0CB00') +theme_void() +
  scale_fill_manual(values= c('lemonchiffon2',  'orangered1',  'chocolate4'), 
                    name= 'Resource Capital Intensity (Ordinal)')+
  theme(legend.title = element_text(size = 12, face = 'bold'),
        legend.text = element_text(size = 14, face = 'bold'))  

mix<- ggarrange(a,b,d,f, labels = c('1830', '1880', '1920', '1950'), 
                ncol=2, nrow=2, legend="right", 
                legend.grob = get_legend(d), common.legend = T)
mix
# ggsave(plot=mix, 'capintsix_color.png', width = 10, height = 6) 




#-------------------------------------------------------
######  Figure 6 ######
#-------------------------------------------------------
 
a=  tapply(base_exp$disputed, base_exp$capord, mean) %>% as.matrix
b=  tapply(uoabord_exp$disputed, uoabord_exp$capord, mean) %>% as.matrix
vec= c('Low', "Med", 'High')
a= cbind(vec,a,b)
colnames(a)= c('bencon', 'claimed', 'highrisk') 
a= as.data.frame(a)
a$claimed <- a$claimed %>% as.numeric()  %>% round(3)
a$claimed_per <- paste0(a$claimed*100, '%')
a$bencon <- factor(a$bencon, levels =    c('Low', "Med", 'High'))
a$highrisk <- a$highrisk %>% as.numeric()  %>% round(3)
a$highrisk_per <- paste0(a$highrisk*100, '%')

p1 = ggplot(data= a, aes(x= bencon, y= claimed)) +
  geom_bar(stat="identity", width = 0.7, fill='grey60')+
  geom_text(aes(y=claimed, label= claimed_per), vjust=2.5,
            color="black", size=4)+ 
  ylab('Percentage of Gridcell-Years Claimed') + 
  xlab('')+
  theme_minimal()  +
  ylim(0, 0.18)+
  scale_x_discrete(labels= c("Low" = "Low \n(Barren;Crops;Forestry)", 'Med'= 'Med \n(Minerals)','High'= 'High \n(Oil & Gas)'))+
  theme(axis.title = element_text(size = 14, face = 'bold'),
        axis.text.x = element_text(size = 12, face = 'bold'))  
p1
 # ggsave('barplot.png',  width=5.4, height=4.9)
 

p2= ggplot(data= a, aes(x= bencon, y= highrisk)) +
  geom_bar(stat="identity", width = 0.7, fill='grey30')+
  geom_text(aes(y=highrisk, label= highrisk_per), vjust=3,
            color="white", size=4)+ 
  ylab('')+
  xlab('')+
  theme_minimal()  +
  ylim(0, 0.18)+
  scale_x_discrete(labels= c("Low" = "Low \n(Barren;Crops;Forestry)", 'Med'= 'Med \n(Minerals)','High'= 'High \n(Oil & Gas)'))+
  theme(axis.title = element_text(size = 14, face = 'bold'),
        axis.text.x = element_text(size = 12, face = 'bold')) 
p2
 # ggsave('barplot_high.png', width=5.4, height=4.9)



#-------------------------------------------------------
###### Table 2 ######
#-------------------------------------------------------

naics <- read_excel("NAICS_capital_primary.xlsx")
naics %<>% filter(`Asset Category`== 'Structures' & `Duration Title`== 'Levels' &
                    `Measure Title`== 'Asset share in industry capital income')
naics %>% filter(NAICS %in% c('111112', '113-115' ,'211', '212')) %>% 
  dplyr::select(NAICS, `NAICS Title`, Mean) %>% arrange(Mean) %>% xtable




#-------------------------------------------------------
######   Table 3 ######
## Warning: This section can take up to 1 hour to run ## 

#-------------------------------------------------------



## make sure data is arranged correctly 
# (spml can turn incorrect estimates if arranged in wrong order)
puoab %<>% arrange(gridid, year) 
pbase %<>% arrange(gridid, year)


sar_time_king <-  spml(disputed ~ capord,  data= puoab, 
                       model = 'within',   effect = 'time',
                       index = c('gridid', 'year'),
                       listw = spatial_w_uoab, lag = T  ,  
                       zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  

sar_all_time_king <- spml(disputed ~ capord+ overlap+ nocol+borderdist,  data= pbase, 
                          model = 'within',   effect = 'time',
                          index = c('gridid', 'year'),
                          listw = spatial_w, lag = T  ,  
                          zero.policy = TRUE, method= 'eigen', spatial.error = 'none') 


summary(sar_time_king)
# capordMed  -0.00335815  0.00068025 -4.9367 7.947e-07 ***
# capordHigh -0.00339809  0.00080654 -4.2132 2.518e-05 ***
summary(sar_all_time_king)
# capordMed  -0.00301363466  0.00041236315  -7.3082 0.0000000000002707 ***
#   capordHigh -0.00385382515  0.00052275975  -7.3721 0.0000000000001680 ***
#   overlap     0.00957575699  0.00020961964  45.6816          < 2.2e-16 ***
#   nocol       0.00951168560  0.00025350697  37.5204          < 2.2e-16 ***
#   borderdist -0.00000201284  0.00000019413 -10.3684          < 2.2e-16 ***


###### marginal effects #####
 
##### all gridcells 
pbasemean= pbase$disputed %>% mean %>% round(3) #9.4%
a=splm::impacts(sar_all_time_king, listw = spatial_w, time=10, 
                set.seed(123)) %>% summary(zstats=T, short=T)
## get raw marginal effects
me = a$res$total%>% round(3)
me 
#  -0.038 -0.049  0.121  0.120  0.000 (for numbers on footnote 18)

## get percentage change 
(me[1]/ pbasemean) %>% round(3)*100 
(me[2]/ pbasemean) %>% round(3)*100 
(me[3]/ pbasemean) %>% round(3)*100 
(me[4]/ pbasemean) %>% round(3)*100 
(a$res$total[5]*100/ pbasemean) %>% round(3)*100 # for border distance scale up 100km

 
#### highrisk gridcells 
puoabmean= puoab$disputed %>% mean %>% round(3) #15.5%
b= splm::impacts(sar_time_king, listw = spatial_w_uoab, time=10,
                 set.seed(123)) %>% summary(zstats=T, short=T)
## get percentage change
me = b$res$total
(me[1]/puoabmean) %>% round(3)*100 
(me[2]/puoabmean) %>% round(3)*100  

 

##########################################################
##################### Appendix ###########################
##########################################################


#-------------------------------------------------------

###### APPENDIX B ######

#-------------------------------------------------------


#------------------------------------------------ 
### B.1: resources by Decade  ######
#------------------------------------------------ 


###### Table B1 ###### 

decade<- str_c('d',seq(1860, 1950, 10),'capint')
decadecap<-  str_c('d',seq(1860, 1950, 10),'capint.Freq')

a= lapply(base, table)

## manually put in this data 
#for decades 1830-50 when there was no oil 
a$d1830capint
a$d1840capint
a$d1850capint


#for all decades after 
a[decade]  %>% as.data.frame() %>% dplyr::select(decadecap) %>% t %>% 
  xtable() 


####### Table B2 ###### 

decade<- str_c('d',seq(1860, 1950, 10),'capint')
decadecap<-  str_c('d',seq(1860, 1950, 10),'capint.Freq')

a= lapply(uoabord, table) 

#for decades 1830-50 when there was no oil 
a$d1830capint
a$d1840capint
a$d1850capint
#for all decades after 
a[decade]  %>% as.data.frame()%>% dplyr::select(decadecap) %>% t %>% 
  xtable() 



#--------------------------------------- 
### B.2: Claimed Cells by Decade  ######
#--------------------------------------- 

######  Table B3  

### Disputed Cells by Decade  
decade<- str_c('d',seq(1830, 1950, 10),'ddisp')
decadecap<- str_c(seq(1830, 1950, 10),'s')

a1<- base %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
colnames(a1) <- c(decadecap)
b1<- uoabord %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
c1<- base %>% filter(uoabord!=1) %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
colnames(b1) <- c(decadecap)
colnames(c1) <- c(decadecap)
rbind(b1,c1,a1) %>% t %>% xtable (digits = 0)


### Initiated Claims by Decade 
decade<- str_c('d',seq(1830, 1950, 10),'init')
decadecap<- str_c(seq(1830, 1950, 10),'s')

a1<- base %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
colnames(a1) <- c(decadecap)
b1<- uoabord %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
c1<- base %>% filter(uoabord!=1) %>%  summarise_at(c( decade), sum, na.rm = TRUE) %>% round(2)
colnames(b1) <- c(decadecap)
colnames(c1) <- c(decadecap)
rbind(b1,c1,a1) %>% t %>% xtable (digits = 0)


#--------------------------------------- 
### B3: Justification for Gridcells as UoA  ######
#--------------------------------------- 

## bring in province data 
provdat <- read_csv("partidos_yr.csv")
provdat = provdat %>% mutate(
  begclaim = as.numeric(begclaim),
  endclaim= as.numeric(endclaim)
) 


# expand
time= seq(1830, 2001, 1)
reprow= rep(rownames(provdat), length(time)) %>% as.numeric()
provdat_exp= provdat[reprow, ]
provdat_exp= provdat_exp %>% arrange(OBJECTID) %>%  group_by(OBJECTID) %>% 
  mutate(dup= 1:n())
#transform to year
provdat_exp= provdat_exp %>% mutate(
  year= 1830+ dup -1
) 
provdat_exp  %<>% arrange(provid, year)

###put in dispute information in to province data 

provdat_exp= provdat_exp %>% group_by(provid) %>% 
  mutate(
    start= ifelse(begclaim== year, 1, 0),
    end= ifelse(endclaim== year, 1, 0),
    check= ifelse(year>=begclaim & year <= endclaim, 1, 0),
    check= replace_na(check, 0))
provdat_exp= provdat_exp %>% group_by(provid, year) %>% 
  mutate(disputed = max(check))  


### filter disputes on claimed area  
### for more explanation see rscript-io-makedata-rep 
### (similar exercise using gridcells)

union <- read_csv("union_part.csv")
landarea <- read_csv("landmass_part.csv")
union$landmass <- landarea$Shape_Area[match(union$provid, landarea$provid)]

claimed= union %>% filter(claim!=0 )
claimed= claimed %>% group_by(claimre, provid) %>% mutate(
  percent= Overlaparea/landmass,
  sumoverlap= sum(percent),
  sumoverlap= ifelse(sumoverlap >= 1 , 1, sumoverlap),
  sumoverlap= ifelse(sumoverlap < 0.05 , 0, sumoverlap)
)

tab= claimed %>% distinct(provid,claimre, sumoverlap)
provdat_exp$disparea <-tab$sumoverlap[match(interaction(provdat_exp$provid, provdat_exp$claimre),
                                           interaction(tab$provid, tab$claimre))]
provdat_exp= provdat_exp %>% mutate(
  disparea= ifelse(is.na(disparea)==T, 0, disparea))
summary(provdat_exp$disparea) 
#   0.000   0.000   0.000   0.134   0.000   1.000 

#####  Figure B3 

## create table for histogram
xx<- seq(0, 1, 0.1)
a<- rep(NA,length(xx))

for (i in 1:length(a)) {
  a[i]<- provdat_exp %>% filter(disputed==1 & disparea>= xx[i] & disparea<xx[i]+0.1) %>%
    pull(provid) %>% unique   %>% length
}
parthist= as.data.frame(cbind(xx, a))
parthist$prop <- parthist$a/ sum(parthist$a)  

#under 10 percent, combine 
parthist[2,]<- parthist[1,]+ parthist[2,]
parthist= parthist[-1,]

## generate figure 
fig<- ggplot(data = parthist, aes(x=reorder(as.factor(xx), -xx), y=a) ) +
  geom_bar(stat='identity') +theme_clean()+
  xlab('Percentage of Province Claimed')+
  ylab('Number of Provinces Claimed') + 
  ylim(0,105)+
  geom_text(aes(label=round(prop,2)), hjust=-0.5, size=3.5)+ coord_flip()
fig
# ggsave( 'provinceper.png', fig, width = 7.2, height = 4.3)




#-------------------------------------------------------

################ APPENDIX D #######################

#-------------------------------------------------------



#--------------------------------------- 
###### D1: decade spatial lag  ######

# *****  WARNING  ****: this table may take up to 8hrs to replicate in full
## please plan accordingly before running whole section

#--------------------------------------- 


## all gridcells 
sp1830cap <- lagsarlm(d1830ddisp ~ d1830capord  +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1840cap <- lagsarlm(d1840ddisp ~ d1840capord  +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1850cap <- lagsarlm(d1850ddisp ~ d1850capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1860cap <- lagsarlm(d1860ddisp ~ d1860capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1870cap <- lagsarlm(d1870ddisp ~ d1870capord +overlap+ nocol+ borderdist, data= base,  
                      listw = spatial_w, zero.policy = T) 
sp1880cap <- lagsarlm(d1880ddisp ~ d1880capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1890cap <- lagsarlm(d1890ddisp ~ d1890capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T) 
sp1900cap <- lagsarlm(d1900ddisp ~ d1900capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)
sp1910cap <- lagsarlm(d1910ddisp ~ d1910capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)
sp1920cap <- lagsarlm(d1920ddisp ~ d1920capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)
sp1930cap <- lagsarlm(d1930ddisp ~ d1930capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)
sp1940cap <- lagsarlm(d1940ddisp ~ d1940capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)
sp1950cap <- lagsarlm(d1950ddisp ~ d1950capord +overlap+ nocol+ borderdist, data= base, 
                      listw = spatial_w, zero.policy = T)


 

### matrix for coefficients 
a1= summary(sp1830cap)$coefficients
a2=   summary(sp1840cap)$coefficients
a3=   summary(sp1850cap)$coefficients
a4=   summary(sp1860cap)$coefficients
a5=   summary(sp1870cap)$coefficients
a6=   summary(sp1880cap)$coefficients
a7=   summary(sp1890cap)$coefficients
a8=   summary(sp1900cap)$coefficients
a9=   summary(sp1910cap)$coefficients
a10=   summary(sp1920cap)$coefficients
a11=   summary(sp1930cap)$coefficients
a12=   summary(sp1940cap)$coefficients
a13=   summary(sp1950cap)$coefficients

decadecap= str_c(seq(1860, 1950, 10), 's')
## for years after 1860s (with oil)
mat = rbind( a4,a5, a6, a7, a8, a9, a10, a11, a12, a13)
row.names(mat) <-decadecap
mat %>% xtable(digits = 3)
## for years before 1860s (no oil)
rbind(a1, a2, a3)%>% xtable(digits = 3)

 

######## For high risk gridcells ##########


sp1830cap_uoa <- lagsarlm(d1830ddisp ~ d1830capord  , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1840cap_uoa <- lagsarlm(d1840ddisp ~ d1840capord  , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1850cap_uoa <- lagsarlm(d1850ddisp ~ d1850capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1860cap_uoa <- lagsarlm(d1860ddisp ~ d1860capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1870cap_uoa <- lagsarlm(d1870ddisp ~ d1870capord , data= uoabord,  
                          listw = spatial_w_uoab, zero.policy = T) 
sp1880cap_uoa <- lagsarlm(d1880ddisp ~ d1880capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1890cap_uoa <- lagsarlm(d1890ddisp ~ d1890capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T) 
sp1900cap_uoa <- lagsarlm(d1900ddisp ~ d1900capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)
sp1910cap_uoa <- lagsarlm(d1910ddisp ~ d1910capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)
sp1920cap_uoa <- lagsarlm(d1920ddisp ~ d1920capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)
sp1930cap_uoa <- lagsarlm(d1930ddisp ~ d1930capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)
sp1940cap_uoa <- lagsarlm(d1940ddisp ~ d1940capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)
sp1950cap_uoa <- lagsarlm(d1950ddisp ~ d1950capord , data= uoabord, 
                          listw = spatial_w_uoab, zero.policy = T)



## matrix for coefficients
a1= summary(sp1830cap_uoa)$coefficients
a2=   summary(sp1840cap_uoa)$coefficients
a3=   summary(sp1850cap_uoa)$coefficients
a4=   summary(sp1860cap_uoa)$coefficients
a5=   summary(sp1870cap_uoa)$coefficients
a6=   summary(sp1880cap_uoa)$coefficients
a7=   summary(sp1890cap_uoa)$coefficients
a8=   summary(sp1900cap_uoa)$coefficients
a9=   summary(sp1910cap_uoa)$coefficients
a10=   summary(sp1920cap_uoa)$coefficients
a11=   summary(sp1930cap_uoa)$coefficients
a12=   summary(sp1940cap_uoa)$coefficients
a13=   summary(sp1950cap_uoa)$coefficients

decadecap= str_c(seq(1860, 1950, 10), 's')
## for years after 1860s (with oil)
mat = rbind( a4,a5, a6, a7, a8, a9, a10, a11, a12, a13)
row.names(mat) <-decadecap
mat %>% xtable(digits = 3)
## for years before 1860s (no oil)
rbind(a1, a2, a3)%>% xtable(digits = 3)
 
  


#--------------------------------
#### D.2: Other Cutoffs ######
# ****Warning*****: this section may take 2-3 hrs to run

#--------------------------------

 
###### robustness checks with 50%  ########

spml(disputed50 ~ capord+ overlap+ nocol+borderdist,  data= pbase, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary

# capordMed  -0.00882511041  0.00041192167 -21.424          < 2.2e-16 ***
#   capordHigh -0.00374487443  0.00052215164  -7.172 0.0000000000007391 ***
#   overlap     0.00486150843  0.00020787922  23.386          < 2.2e-16 ***
#   nocol       0.00803632754  0.00025251523  31.825          < 2.2e-16 ***
#   borderdist -0.00000383826  0.00000019426 -19.758          < 2.2e-16 ***


spml(disputed50 ~ capord,  data= puoab, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w_uoab, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  %>% summary
# capordMed  -0.01445510  0.00068289 -21.1674    < 2.2e-16 ***
#   capordHigh -0.00415422  0.00080956  -5.1315 0.0000002875 ***

spml(disputed50 ~ capord,  data= pbase, 
     model = 'within',   effect = 'twoway',
     index = c('gridid', 'year'),
     listw = spatial_w, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary
# capordMed  -0.00239053  0.00096050 -2.4888 0.012816 * 
#   capordHigh -0.00129532  0.00050082 -2.5864 0.009699 **

spml(disputed50 ~ capord,  data= puoab  , 
     model = 'within',   effect = 'twoway',
     index = c('gridid', 'year'),
     listw = spatial_w_uoab, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary
# capordMed  -0.00020702  0.00160206 -0.1292   0.8972
# capordHigh -0.00013797  0.00078405 -0.1760   0.8603




######## robustness checks with 80% #########


spml(disputed80 ~ capord+ overlap+ nocol+borderdist,  data= pbase, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary
# capordMed  -0.0066241806  0.0004058770 -16.3207   < 2.2e-16 ***
#   capordHigh -0.0023337545  0.0005144807  -4.5361 0.000005729 ***
#   overlap     0.0022828414  0.0002044890  11.1636   < 2.2e-16 ***
#   nocol       0.0064379089  0.0002482279  25.9355   < 2.2e-16 ***
#   borderdist -0.0000043928  0.0000001915 -22.9390   < 2.2e-16 ***


spml(disputed80 ~ capord,  data= puoab, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w_uoab, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  %>% summary
# capordMed  -0.01174752  0.00067069 -17.5157 < 2.2e-16 ***
#   capordHigh -0.00262500  0.00079509  -3.3015 0.0009617 ***

spml(disputed80 ~ capord,  data= pbase, 
     model = 'within',   effect = 'twoway',
     index = c('gridid', 'year'),
     listw = spatial_w, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary
# capordMed  -0.00131143  0.00095152  -1.3782   0.1681    
# capordHigh -0.00528973  0.00049625 -10.6595   <2e-16 ***

spml(disputed80 ~ capord,  data= puoab  , 
     model = 'within',   effect = 'twoway',
     index = c('gridid', 'year'),
     listw = spatial_w_uoab, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary
# capordMed   0.00174754  0.00158166  1.1049             0.2692    
# capordHigh -0.00573365  0.00077422 -7.4057 0.0000000000001305 ***



#--------------------------------
#### D.3: OLS Model ######
#--------------------------------

## check if the data ordering remains the same (should be true)
## important to check if you have run regressions using conleyregs before this section
## running conleyreg version 0.1.7 can mess up the data!!!

all.equal(base_exp$disputed, pbase$disputed %>% as.numeric)
all.equal(uoabord_exp$disputed, puoab$disputed %>% as.numeric)
## should be true
## STOP if false and reload data 

a= felm(disputed~ capord+ overlap + nocol+ borderdist | year| 0| gridid, 
        data= base_exp ) 
b= felm(disputed~ capord|  year| 0| gridid, data=  uoabord_exp ) 

c= conleyreg(disputed~ capord + overlap + nocol+ borderdist  | year, 
             data= base_exp, dist_cutoff = 100,
             lat= 'lat', lon='long', unit= 'gridid', time='year') 
d= conleyreg(disputed~ capord| year, 
             data= uoabord_exp, dist_cutoff = 100,
             lat= 'lat', lon='long', unit= 'gridid', time='year') 
stargazer(a,b,c,d)



# %%%%%  Now we will remove the data we used and reload it  %%%%%%
## !!! conleyreg version 0.1.7 can mess up the original data structure !!!!!
## reloading data is crucial 
 
rm(base_exp)
rm(uoabord_exp)

base_exp <- read_csv("main-gridyear.csv")
base_exp %<>% arrange(gridid, year)
base_exp =base_exp  %>% mutate( 
  capord= factor(capord, ordered = F, levels = c( 'Low','Med',  'High')), 
  capord_orig= factor(capord_orig, ordered = F, levels = c('Low', 'Med', 'High')),
  capord_four= factor(capord_four, ordered = F, levels = c('None','Low', 'Med', 'High')),
  leadcapord= factor(leadcapord, ordered = F, levels = c('Low', 'Med', 'High')),
  leadcapord_noNA= factor(leadcapord_noNA, ordered = F, levels = c('Low', 'Med', 'High'))
)
uoabord_exp <- base_exp %>% filter(uoabord ==1 )


#--------------------------------
#### D.4: Claim Onsets as DV  ######
#--------------------------------

# %%%%%  Attention  %%%%%%
# Make sure to use the NEWLY LOADED DATA #
# Conleyreg version 0.1.7 can mess up data structure) 

## check the data  
all.equal(base_exp$disputed, pbase$disputed %>% as.numeric)
all.equal(uoabord_exp$disputed, puoab$disputed %>% as.numeric)
all.equal(base_exp$init, pbase$init %>% as.numeric)
all.equal(uoabord_exp$init, puoab$init %>% as.numeric)
## should all be true
## STOP if FALSE and reload data 



a= felm(init~ capord+ overlap + nocol+ borderdist | year| 0| gridid, 
        data= base_exp ) 
b= felm(init~ capord|  year| 0| gridid, data= uoabord_exp  )  
c= conleyreg(init~ capord + overlap + nocol+ borderdist  | year, 
             data= base_exp, dist_cutoff = 100,
             lat= 'lat', lon='long', unit= 'gridid', time='year') 
d= conleyreg(init~ capord| year, 
             data= uoabord_exp, dist_cutoff = 100,
             lat= 'lat', lon='long', unit= 'gridid', time='year') 
stargazer(a,b,c,d)



# %%%%%% Now we will remove the data we used and reload it %%%%%%
## !!!!  conleyregs can mess up the original data structure !!!!!
## removing and reloading/ resetting to correct data 

rm(base_exp)
rm(uoabord_exp)

base_exp <- read_csv("main-gridyear.csv")
base_exp %<>% arrange(gridid, year)
base_exp =base_exp  %>% mutate( 
  capord= factor(capord, ordered = F, levels = c( 'Low','Med',  'High')), 
  capord_orig= factor(capord_orig, ordered = F, levels = c('Low', 'Med', 'High')),
  capord_four= factor(capord_four, ordered = F, levels = c('None','Low', 'Med', 'High')),
  leadcapord= factor(leadcapord, ordered = F, levels = c('Low', 'Med', 'High')),
  leadcapord_noNA= factor(leadcapord_noNA, ordered = F, levels = c('Low', 'Med', 'High'))
)
uoabord_exp <- base_exp %>% filter(uoabord ==1 )



#--------------------------------
#### D.5. Two way fixed effects  #######
#--------------------------------
 
sar_two_king <- spml(disputed ~ capord,  data= puoab  , 
                     model = 'within',   effect = 'twoway',
                     index = c('gridid', 'year'),
                     listw = spatial_w_uoab, lag = T  ,  
                     zero.policy = TRUE, method= 'eigen', spatial.error = 'none') 
sar_all_two_king <-  spml(disputed ~ capord,  data= pbase, 
                          model = 'within',   effect = 'twoway',
                          index = c('gridid', 'year'),
                          listw = spatial_w, lag = T  ,  
                          zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  
 


summary(sar_two_king)
# capordMed   0.00246761  0.00163121  1.5127   0.1303
# capordHigh -0.00040670  0.00079831 -0.5095   0.6104 

summary(sar_all_two_king)
# capordMed   0.00131563  0.00098202  1.3397   0.1803  
# capordHigh -0.00128035  0.00051205 -2.5004   0.0124 * 

impacts(sar_two_king, listw = spatial_w_uoab, time=10) %>% summary(zstats=T, short=T)
impacts(sar_all_two_king, listw = spatial_w, time=10) %>% summary(zstats=T, short=T)
 
 

# ########## Table D7 ################

# MAKE SURE YOU ARE USING THE NEWLY LOADED DATA #
# NOT SOMETHING THAT WENT THROUGH CONLEYREG (which can mess up data structure) !!!! 

## check the data (should all be true)
all.equal(base_exp$disputed, pbase$disputed %>% as.numeric)
all.equal(uoabord_exp$disputed, puoab$disputed %>% as.numeric)
all.equal(base_exp$init, pbase$init %>% as.numeric)
all.equal(uoabord_exp$init, puoab$init %>% as.numeric)
## should all be true
## STOP if FALSE and reload data 




a= felm(disputed~ capord| year +gridid | 0| gridid,
        data= base_exp )
b= felm(disputed~ capord|  year +gridid | 0| gridid, data= uoabord_exp  )
c= felm(init~ capord | year +gridid | 0| gridid,
        data= base_exp )
d= felm(init~ capord|  year +gridid | 0| gridid, data= uoabord_exp  )

stargazer(a,b,c,d)
  

#----------------------------------------------------------
#### D.6. Using only original USGS mineral data #######
#----------------------------------------------------------


sarorig<- spml(disputed ~ capord_orig,  data= puoab, 
               model = 'within',   effect = 'time',
               index = c('gridid', 'year'),
               listw = spatial_w_uoab, lag = T  ,  
               zero.policy = TRUE, method= 'eigen', spatial.error = 'none') 
summary(sarorig)
# capord_origMed  -0.00843290  0.00109413 -7.7074 1.284e-14 ***
# capord_origHigh -0.00342697  0.00080619 -4.2508 2.130e-05 ***
 
sarorig_two <- spml(disputed ~ capord_orig,  data= puoab  , 
                    model = 'within',   effect = 'twoway',
                    index = c('gridid', 'year'),
                    listw = spatial_w_uoab, lag = T  ,  
                    zero.policy = TRUE, method= 'eigen', spatial.error = 'none') 
summary(sarorig_two)
# capord_origMed  -0.00418029  0.00257517 -1.6233   0.1045
# capord_origHigh -0.00072285  0.00078222 -0.9241   0.3554

sarorig_all <- spml(disputed ~ capord_orig+ overlap+ nocol+borderdist,  data= pbase, 
                    model = 'within',   effect = 'time',
                    index = c('gridid', 'year'),
                    listw = spatial_w, lag = T  ,  
                    zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  
summary(sarorig_all)
# capord_origMed  -0.00558023239  0.00069634169  -8.0136 1.114e-15 ***
# capord_origHigh -0.00382766074  0.00052259351  -7.3244 2.400e-13 ***
#   overlap          0.00958644577  0.00020963227  45.7298 < 2.2e-16 ***
#   nocol            0.00952121366  0.00025345701  37.5654 < 2.2e-16 ***
#   borderdist      -0.00000199261  0.00000019411 -10.2653 < 2.2e-16 ***

sarorig_two_all <- spml(disputed ~ capord,  data= pbase,
                        model = 'within',   effect = 'twoway',
                        index = c('gridid', 'year'),
                        listw = spatial_w, lag = T  ,
                        zero.policy = TRUE, method= 'eigen', spatial.error = 'none')
summary(sarorig_two_all)
# capordMed   0.00131563  0.00098202  1.3397   0.1803  
# capordHigh -0.00128035  0.00051205 -2.5004   0.0124 *


## calculate marginal effects
impacts(sarorig, listw = spatial_w_uoab, time=10) %>% summary(zstats=T, short=T)
impacts(sarorig_all, listw = spatial_w, time=10) %>% summary(zstats=T, short=T)
impacts(sarorig_two, listw = spatial_w_uoab, time=10) %>% summary(zstats=T, short=T)
impacts(sarorig_two_all, listw = spatial_w, time=10) %>% summary(zstats=T, short=T)



#--------------------------------
#### D.7. Hazard Model #######
#--------------------------------
 
fit <- coxph(Surv(year, I(year+1), disputed ) ~ capord + overlap+ nocol+ borderdist+
               cluster(gridid), data=base_exp)
summary(fit)




#--------------------------------
####### D8. Terrain ##########
#--------------------------------

## read in terrain data 
terrain <- read_excel("terrain.xlsx")

## make neighboring matrix for data with terrain information 

nb<- poly2nb(visual %>%  filter(OBJECTID  %in% terrain$Grid ) , queen = F)
spatial_w_terrain <- nb2listw(nb, style = 'W', zero.policy = T)  

nb<- poly2nb(uoabord_vis %>%  filter(OBJECTID  %in% terrain$Grid ) , queen = F)
spatial_w_uoa_terrain <- nb2listw(nb, style = 'W', zero.policy = T)  


## spatial regressions 

spml(disputed ~ capord + log(terr_mean+1) + log(terr_range+1), 
     data= puoab %>% filter(gridid  %in% terrain$Grid ), 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w_uoa_terrain, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none') %>% summary

#   lambda 0.92165536 0.00043913  2098.8 < 2.2e-16 ***
#   Estimate  Std. Error t-value  Pr(>|t|)    
# capordMed           -0.00149131  0.00069585 -2.1431           0.032102 *  
#   capordHigh          -0.00239046  0.00081149 -2.9458           0.003221 ** 
#   log(terr_mean + 1)  -0.00024640  0.00018146 -1.3579           0.174504    
# log(terr_range + 1) -0.00130500  0.00017802 -7.3307 0.0000000000002289 ***


# 
spml(disputed ~ capord+ overlap+ nocol+borderdist+ log(terr_mean+1) +log(terr_range+1),  
     data= pbase %>% filter(gridid  %in% terrain$Grid ), 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w_terrain, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none') %>% summary

# lambda 0.93019522 0.00031609  2942.8 < 2.2e-16 ***
# 
# Coefficients:
# capordMed           -0.00189449142  0.00042008950 -4.5097 6.491e-06 ***
#   capordHigh          -0.00332679417  0.00052689599 -6.3139 2.720e-10 ***
#   overlap              0.00976580973  0.00021172270 46.1255 < 2.2e-16 ***
#   nocol                0.00918138862  0.00025462106 36.0590 < 2.2e-16 ***
#   borderdist          -0.00000157291  0.00000019677 -7.9937 1.310e-15 ***
#   log(terr_mean + 1)  -0.00008465359  0.00010629379 -0.7964    0.4258    
# log(terr_range + 1) -0.00097312844  0.00011417715 -8.5230 < 2.2e-16 ***





#-------------------------------------------------------

###### APPENDIX E ######

#-------------------------------------------------------
 


#----------------------------------------------------------
###### E.1: Previously settled disputes #########
#----------------------------------------------------------


sar_time_settled <-  spml(disputed ~ capord+settled,  data= puoab, 
                          model = 'within',   effect = 'time',
                          index = c('gridid', 'year'),
                          listw = spatial_w_uoab, lag = T  ,  
                          zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  
summary(sar_time_settled)

#   Estimate  Std. Error  t-value     Pr(>|t|)    
# capordMed  -0.00357305  0.00068023  -5.2527 0.0000001499 ***
#   capordHigh -0.00345965  0.00080644  -4.2900 0.0000178643 ***
#   settled    -0.00785228  0.00038350 -20.4755    < 2.2e-16 ***



sar_all_time_settled <-  spml(disputed ~ capord+ overlap+ nocol+borderdist+settled,  data= pbase, 
                              model = 'within',   effect = 'time',
                              index = c('gridid', 'year'),
                              listw = spatial_w, lag = T  ,  
                              zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  
summary(sar_all_time_settled)

#   Estimate     Std. Error  t-value  Pr(>|t|)    
# capordMed  -0.00307640509  0.00041228158  -7.4619 8.528e-14 ***
#   capordHigh -0.00390279928  0.00052265650  -7.4672 8.190e-14 ***
#   overlap     0.01016734650  0.00021063889  48.2691 < 2.2e-16 ***
#   nocol       0.01189589569  0.00026543434  44.8167 < 2.2e-16 ***
#   borderdist -0.00000239721  0.00000019455 -12.3218 < 2.2e-16 ***
#   settled    -0.00846454736  0.00027272744 -31.0367 < 2.2e-16 ***


## marginal effects 
impacts(sar_time_settled, listw = spatial_w_uoab, time=10) %>% summary(zstats=T, short=T)
impacts(sar_all_time_settled, listw = spatial_w, time=10) %>% summary(zstats=T, short=T)




#----------------------------------------------------------
###### E.2: Anticipated Discoveries #########
#----------------------------------------------------------

spml(disputed ~ leadcapord_noNA,  data= puoab, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w_uoab, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')  %>% summary
# lambda 0.92208902 0.00043761  2107.1 < 2.2e-16 ***
# leadcapord_noNAMed  -0.00379604  0.00068126 -5.5721 0.00000002517 ***
#   leadcapord_noNAHigh -0.00318480  0.00074474 -4.2764 0.00001899481 ***

spml(disputed ~ leadcapord_noNA+ overlap+ nocol+borderdist,  data= pbase, 
     model = 'within',   effect = 'time',
     index = c('gridid', 'year'),
     listw = spatial_w, lag = T  ,  
     zero.policy = TRUE, method= 'eigen', spatial.error = 'none')   %>% summary 
# lambda 0.92971903 0.00031331  2967.4 < 2.2e-16 ***
# 
# Coefficients:
#   Estimate     Std. Error  t-value  Pr(>|t|)    
# leadcapord_noNAMed  -0.00318509842  0.00041230107  -7.7252 1.117e-14 ***
#   leadcapord_noNAHigh -0.00346585456  0.00047966124  -7.2256 4.988e-13 ***
#   overlap              0.00958183616  0.00020964342  45.7054 < 2.2e-16 ***
#   nocol                0.00950784483  0.00025356410  37.4968 < 2.2e-16 ***
#   borderdist          -0.00000201162  0.00000019414 -10.3615 < 2.2e-16 ***
 
 
 
########### Table E3 ################


a= felm(resourcedisc ~ settled|year|0|gridid, pbase) 
b= felm(resourcedisc ~ settled|year|0|gridid, pbase %>%
          filter(name != 'Patagonia'| is.na((name)==T))) 
stargazer(a, b)



#----------------------------------------------

############### Section E.3 ###################

#----------------------------------------------

######## Preparation for Section E.3 #############

### for this section we need to call in dyad-year-grid data
fin<- read_csv("dyad-year-gridcell.csv")
fin$capord = factor(fin$capord, levels = c('Low', 'Med', 'High') )


#----------------------------------------------------------------------
####  run manual interaction plot function 
interfelm = function(reg1=reg1, x1= x1, x2= x2, data= data){
  beta.hat <- coef(reg1)
  vcov1 <- vcov(reg1)
  interterm= paste0(x1,':',x2)
  datax2 = data [ ,x2]
  datax1 = data [ ,x1]
  z0 <- seq(min(datax2, na.rm = T), max(datax2, na.rm = T),
            length.out = 100)
  dy.dx <- beta.hat[x1] + beta.hat[interterm]*z0
  se.dy.dx <- sqrt(vcov1[x1, x1] + (z0^2)*vcov1[interterm, interterm] + 
                     2*z0*vcov1[x1, interterm])
  upr <- dy.dx + 1.96*se.dy.dx
  lwr <- dy.dx - 1.96*se.dy.dx
  
  a=  ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
    labs(x=x2, y="Marginal Effects",
         # title=paste("Marginal Effects of", x1, "on Y by",  x2), 
         cex=4) +
    geom_line(aes(z0, dy.dx),lwd = 1) +
    geom_hline(yintercept=0, lwd = 1, linetype=3) +
    geom_ribbon(aes(ymin=lwr, ymax=upr), alpha=0.3)
  return(a)
}

 
#-----------------------------------------------------------------


########## Table E4 ############

a= felm(realdisp~ capord  | year +dyadid | 0| gridid, data= fin %>% filter(bopimbal>=0.25) )  
b= felm(realinit~ capord | year +dyadid | 0| gridid, data= fin %>% filter(bopimbal>=0.25) )   
c= felm(realdisp~ capord  | year  | 0| gridid, data= fin %>% filter(bopimbal>=0.25) )  
d= felm(realinit~ capord | year  | 0| gridid, data= fin %>% filter(bopimbal>=0.25) )   
stargazer(c,d, a,b)


########## Table E5 #############
a= felm(realdisp~ capord  | year +dyadid | 0| gridid, data= fin %>% filter( colonial ==0 &   bopimbal>=0.25) )  
b= felm(realinit~ capord | year +dyadid | 0| gridid, data= fin %>% filter( colonial ==0 & bopimbal>=0.25) )   
c= felm(realdisp~ capord  | year  | 0| gridid, data= fin %>% filter( colonial ==0 & bopimbal>=0.25) )  
d= felm(realinit~ capord | year  | 0| gridid, data= fin %>% filter( colonial ==0 & bopimbal>=0.25) )   
stargazer(c,d, a,b)



######### Figure E1 ##########

#### E1(a)
reg1 = felm(realdisp~ capord_bin*bopimbal | year+dyadid| 0| gridid, data= fin )
interfelm(reg1, x1= 'capord_bin', x2='bopimbal', fin)   +   
  geom_density(data=fin, aes(x= bopimbal, y= ..scaled../30 -0.08), inherit.aes = F, alpha = 0.4) +
  labs(title ='Effect of resource presence by Balance of Power') +xlab('Power Imbalance' )
# ggsave('bopdisp.png',width = 7.5, height = 5.4)


#### E1(b)
reg1 = felm(realinit~ capord_bin*bopimbal | year +dyadid | 0| gridid, data= fin )
interfelm(reg1, x1= 'capord_bin', x2='bopimbal', fin)    +   
  geom_density(data=fin, aes(x= bopimbal, y= ..scaled../1000 -0.002 ), inherit.aes = F, alpha = 0.4) +
  labs(title ='Effect of resource presence by Balance of Power') +xlab('Power Imbalance' )
# ggsave('bopinit.png', width = 7.5, height = 5.4)




########## Table E6 ############
a= felm(realdisp~ capord | year +dyadid | 0| gridid, data= fin %>% filter(tenchange ==1) )   
b= felm(realinit~ capord | year +dyadid | 0| gridid, data= fin %>% filter(tenchange ==1) )  
c= felm(realdisp~ capord | year +dyadid | 0| gridid, data= fin %>% filter(twentychange ==1) )   
d= felm(realinit~ capord | year +dyadid | 0| gridid, data= fin %>% filter(twentychange ==1) )   
stargazer(a,b,c,d)
 



 
 #--------------------------------------------------------
 
 ######### Appendix G ############## 
 
 #--------------------------------------------------------

 

 
 ########### Table G1 #############
 
 #### Call in data 
 grids<- read_csv('dyad-grid.csv')
 grids %<>% mutate(
   dyname = case_when(dyadid== 100101 ~ 'COL-VEN',
                      dyadid== 100130 ~ 'COL-ECU',
                      dyadid== 100135 ~ 'COL-PER',
                      dyadid== 100140 ~ 'COL-BRA',
                      dyadid== 101110 ~ 'VEN-GUY',
                      dyadid== 101140 ~ 'VEN-BRA',
                      dyadid== 110115 ~ 'GUY-SUR',
                      dyadid== 110140 ~ 'GUY-BRA',
                      dyadid== 115140 ~ 'SUR-BRA',
                      dyadid== 115220 ~ 'SUR-FRG',
                      dyadid== 130135 ~ 'ECU-PER',
                      dyadid== 135140 ~ 'PER-BRA',
                      dyadid== 135145 ~ 'PER-BOL',
                      dyadid== 135155 ~ 'PER-CHL',
                      dyadid== 140145 ~ 'BRA-BOL',
                      dyadid== 140150 ~ 'BRA-PAR',
                      dyadid== 140160 ~ 'BRA-ARG',
                      dyadid== 140165 ~ 'BRA-URG',
                      dyadid== 140220 ~ 'BRA-FRG',
                      dyadid== 145150 ~ 'BOL-PAR',
                      dyadid== 145155 ~ 'BOL-CHL',
                      dyadid== 145160 ~ 'BOL-ARG',
                      dyadid== 150160 ~ 'PAR-ARG',
                      dyadid== 155160 ~ 'CHL-ARG',
                      dyadid== 160165 ~ 'ARG-URG',
                      T~ NA
   )
 ) 
 tapply(grids$gridid, grids$dyname, function(x) length(unique(x))) %>% 
   data.frame()%>%  xtable() 
 
 
 
 
 ############ Figure G1 ############## 
 
 vdem_SA<- read_csv('vdem_SA.csv')
 
 vdem_SA %>% filter(year>=1830 & year <2002)  %>% ggplot(aes(year, y= v2psoppaut_osp))+
   geom_line()+ facet_wrap(~country_name, ncol=3)+ theme(axis.title =  element_blank())+
   ylim(0,4)
 # ggsave('oppaut-time.png', width= 9.2, height = 8.2)
 
 
 
 ############ Figure G2 ############## 
 
 
 vdem_SA %>% filter(year>=1830 & year <2002) %>% ggplot(aes(year, y= v2pepwrses_osp))+
   geom_line()+ facet_wrap(~country_name, ncol=3)+ theme(axis.title =  element_blank())+
   ylim(0,4)
 # ggsave('soceq-time.png', width= 9.2, height = 8.2)
 
 
 ########### Table G2 ################
 
 a= felm(realdisp~ capord_bin*oppaut_av | year+ dyadid | 0| dyadid, data= fin )    
 b= felm(realdisp~ capord_bin*soceq_av  | year+ dyadid | 0| dyadid, data= fin )   
 c= felm(realdisp~ capord_bin*oppaut_av+  capord_bin*soceq_av +bop +prevclaimed | year+ dyadid | 0| dyadid, data= fin ) 
 d= felm(realinit~ capord_bin*oppaut_av | year+ dyadid | 0| dyadid, data= fin )    
 e= felm(realinit~ capord_bin*soceq_av  | year+ dyadid | 0| dyadid, data= fin )   
 f= felm(realinit~ capord_bin*oppaut_av+  capord_bin*soceq_av +bop +prevclaimed | year+ dyadid | 0| dyadid, data= fin ) 
 stargazer(a,b,c,d, e,f)
 
 
 
 
 ############ Figure G3 ##############
 
 reg1 =  felm(realdisp~ capord_bin*oppaut_av | year  | 0| dyadid, data= fin ) 
 interfelm(reg1, x1= 'capord_bin', x2='oppaut_av', fin)  +
   geom_density(data=fin, aes(x= oppaut_av, y= ..scaled../15- 0.18), 
                inherit.aes = F, alpha = 0.4, adjust=2.5)+
   labs(title ='Effect of Resource Presence by Opposition Strength') +
   xlab('Opposition Autonomy (Average Dyad Value)')+
   theme(axis.title  = element_text(size=17), title  = element_text(size=17),
         axis.text = element_text(size=12))
  # ggsave('oppdisp-av.png', width = 8.5, height = 5.5)
 
 reg1 = felm(realinit~ capord_bin*oppaut_av | year | 0| dyadid, data= fin )
 interfelm(reg1, x1= 'capord_bin', x2='oppaut_av', fin)  +   
   geom_density(data=fin, aes(x= oppaut_av, y= ..scaled../900- 0.004), 
                inherit.aes = F, alpha = 0.4, adjust= 2.5)+ 
   labs(title ='Effect of Resource Presence by Opposition Strength') +
   xlab('Opposition Autonomy (Average Dyad Value)')+
   theme(axis.title  = element_text(size=17), title  = element_text(size=17),
         axis.text = element_text(size=12))
 # ggsave('oppinit-av.png', width = 8.5, height = 5.5)
 
 
 
 ############ Figure G4 ##############
 
 
 reg1 = felm(realdisp~ capord_bin*soceq_av | year | 0| dyadid, data= fin )
 interfelm(reg1, x1= 'capord_bin', x2='soceq_av', fin)  +
   geom_density(data=fin, aes(x= soceq_av, y= ..scaled../20- 0.15), 
                inherit.aes = F, alpha = 0.4, adjust=2)+
   labs(title ='Effect of Resource Presence by Credibility of Redistribution') +
   xlab('Socioeconomic Equality (Average Dyad Value)')+
   theme(axis.title  = element_text(size=17), title  = element_text(size=17),
         axis.text = element_text(size=12))
  # ggsave('soceqdisp-av.png', width = 8.5, height = 5.5)
 
 
 reg1 = felm(realinit~ capord_bin*soceq_av | year | 0| dyadid, data= fin )
 interfelm(reg1, x1= 'capord_bin', x2='soceq_av', fin)  +   
   geom_density(data=fin, aes(x= soceq_av, y= ..scaled../700- 0.003),
                inherit.aes = F, alpha = 0.4, adjust= 2) +
   labs(title ='Effect of Resource Presence by Credibility of Redistribution') +
   xlab('Socioeconomic Equality (Average Dyad Value)')+
   theme(axis.title  = element_text(size=17), title  = element_text(size=17),
         axis.text = element_text(size=12))
  # ggsave('soceqinit-av.png', width = 8.5, height = 5.5)
 
 
 
 #------------------------------------------------------------------------------
 
 ############################### END of script ############################### 
 
 #-----------------------------------------------------------------------------
 
 
